This file evaluates seagrass gains and losses in coverage from 2018 to 2020. The goal is to characterize the distribution of depth estimates for areas where changes have occurred. A hypothesis is that seagrass may be lost at deeper depths as compared to where it was gained, possibly related to changes in water clarity that can cause seagrass to grow in more shallow areas where they are less light limited. Source content is here.

The analysis is in several steps:

  1. A change analysis between 2018 and 2020 to identify areas where seagrass occurred in 2018 but not in 2020 (loss) and areas where seagrass did not occur in 2018 but did occur in 2020 (gain).
  2. A bathymetry layer is then used to estimate the average depth, as meters below MLLW, for each polygon identified as lost or gained.
  3. The summarized polygons are then intersected with bay segment and seagrass management areas to evaluate potential changes by major areas in the bay.

First, the relevant libraries and datasets are imported. All analyses are conducted using the NAD83(HARN) / Florida West (ftUS) projection.

knitr::opts_chunk$set(echo = TRUE, warning = F, message = F)

library(sf)
library(tidyverse)
library(tbeptools)
library(mapview)
library(units)
library(stars)
library(raster)
library(here)

data(dem)
data(sgseg)
data(sgdat2018)
data(sgdat2020)

# this is the projstring for NAD83(HARN) / Florida West (ftUS)
# https://spatialreference.org/ref/epsg/2882/
prj4 <- '+proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs'

# reproject layers
tbseg <- tbseg %>% 
  st_transform(crs = prj4)
sgmanagement <- sgmanagement %>% 
  st_transform(crs = prj4)

# reproject dem, convert to stars 
demstr <- dem %>% 
  projectRaster(crs = prj4) %>% 
  st_as_stars()

# reproject seagrass layers
sgdat2018 <- sgdat2018 %>% 
  st_transform(crs = prj4)
stdat2020 <- sgdat2020 %>% 
  st_transform(crs = prj4)

This code chunk does the following:

  1. Filters continuous and patchy seagrass from the 2018 and 2020 layers, unions by features
  2. Conducts a true change analysis by taking the differences between the 2018 and 2020 seagrass layers, then their intersection
  3. Combines all results showing changes between seagrass categories (source as 2018 category, target as 2020 category)
  4. Saves the results because the change analysis takes a while.
# 2018 seagrass data
a <- sgdat2018 %>% 
  left_join(fluccs, by = 'FLUCCSCODE') %>%
  st_union(by_feature = TRUE) %>%
  mutate(
    Category = paste0(Category, ', 2018')
  )

# 2020 seagrass data
b <- sgdat2020 %>% 
  st_union(by_feature = TRUE) %>%
  left_join(fluccs,by = 'FLUCCSCODE') %>% 
  mutate(
    Category = paste0(Category, ', 2020')
  )
  
# so intersect doesnt complain about attributes
st_agr(a) = "constant"
st_agr(b) = "constant"

# union separate layers for faster intersect
aunion <- a %>% 
  st_union %>% 
  st_set_precision(1e5) %>% 
  st_make_valid() %>% 
  st_buffer(dist = 0)
bunion <- b %>% 
  st_union %>% 
  st_set_precision(1e5) %>% 
  st_make_valid() %>% 
  st_buffer(dist = 0)
  
# get full union
op1 <- st_difference(a, bunion)
op2 <- st_difference(b, aunion) %>%
  rename(Category.1 = Category)
op3 <- st_intersection(a, b)

# combine and make multipolygon to polygon, do by id otherwise it doesn't work
unidat <- bind_rows(op1, op2, op3) %>%
  mutate(
    yr = unique(na.omit(sub('^.*\\,\\s([0-9]+)$', '\\1', Category))),
    yr.1 = unique(na.omit(sub('^.*\\,\\s([0-9]+)$', '\\1', Category.1))),
    Category = ifelse(is.na(Category), paste0('other, ', yr), as.character(Category)),
    Category.1 = ifelse(is.na(Category.1), paste0('other, ', yr.1), as.character(Category.1)),
    idval = 1:nrow(.)
  ) %>%
  dplyr::select(idval, Category.1, Category) %>%
  dplyr::select(idval, source = Category, target = Category.1) %>% 
  group_by(idval, source, target) %>% 
  nest() %>% 
  mutate(
    data = purrr::map(data, st_cast, 'POLYGON')
  ) %>% 
  unnest('data') %>% 
  st_as_sf()

The unidat dataset is subset to only gains and losses (i.e., polygons can remain the same between years or change between seagrass categories). Areas of each polygon are calculated in acres. Polygons less than 0.023 acres are also removed because these are smaller than the cell size (30 by 30 feet, or 0.023 acres) of the bathymetry layer. The file is saved because it takes a long time to calculate.

# select only lost/gained
chgdat <- unidat %>% 
  ungroup() %>% 
  filter(target == 'other, 2020' | source == 'other, 2018') %>% 
  mutate(
    Acres = st_area(.), 
    Acres = set_units(Acres, 'acres'), 
    Acres = as.numeric(Acres), 
    var = case_when(
      target == 'other, 2020' ~ 'lost', 
      source == 'other, 2018' ~ 'gained'
    )
  ) %>% 
  dplyr::filter(Acres > 0.023) # remove slivers and those less than the pixel size (dem pixel size is about 30x30ft)

save(chgdat, file = here('data/chgdat.RData'))

Reload the data for the rendered file. The results show polygon types in 2018 and 2020 in the source and target columns. A polygon type of ‘other’ in the source column means seagrass was gained in 2020 and a polygon type of ‘other’ in the target column means seagrass was lost in 2020.

load(file = here('data/chgdat.RData'))
chgdat
## Simple feature collection with 4310 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 381530.6 ymin: 1150157 xmax: 531285.6 ymax: 1345319
## CRS:           +proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs
## # A tibble: 4,310 x 6
##    idval source       target                               geometry  Acres var  
##  * <int> <chr>        <chr>              <POLYGON [US_survey_foot]>  <dbl> <chr>
##  1     7 Patchy, 2018 other, 2020 ((411859.5 1250585, 411851.4 125~  0.841 lost 
##  2     8 Patchy, 2018 other, 2020 ((410664.3 1264177, 410671.8 126~  0.237 lost 
##  3     8 Patchy, 2018 other, 2020 ((411085.6 1262250, 411087.4 126~  0.973 lost 
##  4     8 Patchy, 2018 other, 2020 ((411220.3 1261246, 411225.4 126~  0.239 lost 
##  5    10 Patchy, 2018 other, 2020 ((391000.9 1273733, 390998.9 127~  0.869 lost 
##  6    15 Patchy, 2018 other, 2020 ((416303.4 1211170, 416301.8 121~ 15.0   lost 
##  7    16 Patchy, 2018 other, 2020 ((416374.9 1210817, 416377.2 121~  0.248 lost 
##  8    17 Patchy, 2018 other, 2020 ((415834.4 1209104, 415845.7 120~  0.902 lost 
##  9    22 Patchy, 2018 other, 2020 ((385843.3 1274880, 385809.3 127~  0.218 lost 
## 10    25 Patchy, 2018 other, 2020 ((382545.8 1284993, 382537.6 128~  0.109 lost 
## # ... with 4,300 more rows
m <- mapview(dem, layer.name = 'Depth (m)') + mapview(chgdat, zcol = 'var', layer.name = 'Seagrass', col.regions = c("#F8766D", "#00BFC4"))
m


Next, average depth for each polygon is esimated from the bathymetry layer. The chgdat layer (gained/lost) is aggregated by taking the average of the cells from the bathymetry layer that are within each polygon. Polygons are removed that have no mean depth estimate because they are on the boundary or outside of the bathymetry layer. This also takes a while and is not rendered in this file.

chgdatarea <- chgdat %>% 
  aggregate(demstr, ., mean, as_points = T) %>% 
  st_as_sf() %>% 
  rename(meandepth = GDAL.Band.Number.1) %>% 
  bind_cols(st_set_geometry(chgdat, NULL)) %>% 
  filter(!is.na(meandepth))

save(chgdatarea, file = here('data/chgdatarea.RData'))

Reload the data. The meandepth column shows the average depth within each polygon as meters below MLLW (negative values).

load(file = here('data/chgdatarea.RData'))
chgdatarea
## Simple feature collection with 3236 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 409944.4 ymin: 1150701 xmax: 531285.6 ymax: 1345319
## CRS:           +proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs
## First 10 features:
##      meandepth idval       source      target      Acres  var
## 1  -0.69418468    27  Cont., 2018 other, 2020 0.10473300 lost
## 2   0.43311293    31  Cont., 2018 other, 2020 0.04347834 lost
## 3  -1.19202473    49 Patchy, 2018 other, 2020 0.30698955 lost
## 4   0.21248606    49 Patchy, 2018 other, 2020 0.37404466 lost
## 5   0.08625199    49 Patchy, 2018 other, 2020 0.08304300 lost
## 6  -1.04209483    52 Patchy, 2018 other, 2020 0.15276761 lost
## 7  -1.14997934    52 Patchy, 2018 other, 2020 0.21398557 lost
## 8  -1.14677923    52 Patchy, 2018 other, 2020 0.37341319 lost
## 9  -1.26787841    53 Patchy, 2018 other, 2020 2.86281831 lost
## 10 -1.59664105    57 Patchy, 2018 other, 2020 2.40369859 lost
##                          geometry
## 1  POLYGON ((420115.7 1214353,...
## 2  POLYGON ((429474.7 1229396,...
## 3  POLYGON ((439513.2 1225601,...
## 4  POLYGON ((441174.7 1225210,...
## 5  POLYGON ((441836.9 1224821,...
## 6  POLYGON ((426139.4 1204520,...
## 7  POLYGON ((426142.6 1204596,...
## 8  POLYGON ((426455.1 1204737,...
## 9  POLYGON ((426557.5 1204604,...
## 10 POLYGON ((421091.5 1196938,...

Now the change layer with depth estimates is combined with segment and management polygons to summarize by areas of interest.

# intersect change with relevant boundaries
chgdatarea <- chgdatarea %>% 
  st_intersection(tbseg) %>% 
  st_intersection(sgmanagement) %>% 
  unite('allbnds', bay_segment, areas, sep = ', ', remove = F)
chgdatarea
## Simple feature collection with 2620 features and 10 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 416881.4 ymin: 1157848 xmax: 528821.4 ymax: 1345074
## CRS:           +proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs
## First 10 features:
##        meandepth idval       source       target       Acres    var
## 114   0.11488018   629 Patchy, 2018  other, 2020 50.90067840   lost
## 235   0.04889106  1055 Patchy, 2018  other, 2020  9.64494620   lost
## 240  -0.17614305  1059 Patchy, 2018  other, 2020 96.57784211   lost
## 1287 -0.39206555  2729 Patchy, 2018  other, 2020  2.87894694   lost
## 1288 -0.43271106  2729 Patchy, 2018  other, 2020  0.24565248   lost
## 1289 -0.45060265  2729 Patchy, 2018  other, 2020  0.06399648   lost
## 1290 -0.40984065  2729 Patchy, 2018  other, 2020 13.52472502   lost
## 1291 -0.20258051  2729 Patchy, 2018  other, 2020  1.27926961   lost
## 1292 -0.11987968  2729 Patchy, 2018  other, 2020  1.18314739   lost
## 1836 -0.43993417  4222  other, 2018 Patchy, 2020  0.12147729 gained
##             long_name allbnds bay_segment areas                       geometry
## 114  Hillsborough Bay   HB, 2          HB     2 MULTIPOLYGON (((497592.4 13...
## 235  Hillsborough Bay   HB, 2          HB     2 POLYGON ((499438.1 1293914,...
## 240  Hillsborough Bay   HB, 2          HB     2 POLYGON ((503023.2 1306953,...
## 1287 Hillsborough Bay   HB, 3          HB     3 POLYGON ((521900.9 1295906,...
## 1288 Hillsborough Bay   HB, 3          HB     3 POLYGON ((522238.2 1295609,...
## 1289 Hillsborough Bay   HB, 3          HB     3 POLYGON ((522404.3 1294798,...
## 1290 Hillsborough Bay   HB, 3          HB     3 POLYGON ((525674.3 1281558,...
## 1291 Hillsborough Bay   HB, 3          HB     3 POLYGON ((526041.1 1281997,...
## 1292 Hillsborough Bay   HB, 3          HB     3 POLYGON ((526250.7 1283956,...
## 1836 Hillsborough Bay   HB, 3          HB     3 POLYGON ((522074.6 1295830,...

The plot below shows the cumulative distribution of acres lost and gained by major bay segments.

ggplot(chgdatarea, aes(x = meandepth, color = var)) +
  stat_ecdf() +
  facet_wrap(~bay_segment, ncol = 4) +
  theme_minimal() +
  labs(
    color = NULL,
    x = 'Mean depth below MLLW (m)', 
    y = 'ECDF'
  )

mapview(tbseg, legend = F, layer.name = 'Bay segments')


The plot below shows the cumulative distribution of acres lost and gained by major bay segments and seagrass management areas.

ggplot(chgdatarea, aes(x = meandepth, color = var)) +
  stat_ecdf() +
  facet_wrap(~allbnds) +
  theme_minimal() +
  labs(
    color = NULL,
    x = 'Mean depth below MLLW (m)', 
    y = 'ECDF'
  )

mapview(sgmanagement, legend = F, layer.name = 'Management areas', zcol = NULL)